# APPENDIX ---------------------------------------------------------------
# Logit ---------------------------------------------------------------

r2_choices_num_with_lasso <- r2_choices_num %>%
  mutate(
    across(any_of(control_vars),
           ~ .x * as.numeric(pair_includes_trans),
           .names = "{.col}__pair_includes_trans")
  ) %>%
  group_by(group_id) %>%
  mutate(
    across(any_of(control_vars),
           list(group_control = ~ (sum_na(.x) - .x) / (n() - 1)),
           .names = "{.col}__group_control")
  ) %>%
  ungroup()

logits <- list(
  "Chose worker in private\noutcome round (=1)" = glm(
    r2_choose_comparator ~ pair_includes_trans + pair_includes_trans_discussion_full + discussion_full,
    data = r2_choices_num_with_lasso %>% filter(discuss_type %in% c("control", "discussion_full")),
    family = quasibinomial(link = 'logit')
  ),
  "Chose worker in private\noutcome round (=1)" = glm(
    as.formula(
      str_glue("r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_discussion_full + discussion_full + pair_includes_trans_alt:(stratum_id + video_type_placebo + video_type_treatment + delivery_incentive_exp + {model_list[[2]]$lasso_controls %>% fml_sum()}) + (stratum_id + video_type_placebo + video_type_treatment + delivery_incentive_exp + {model_list[[2]]$lasso_controls %>% fml_sum()}) + item_diff + r2_reliability_diff * r2_reliability_shown ")
    ),
    data = r2_choices_num_with_lasso %>% filter(discuss_type %in% c("control", "discussion_full")),
    family = quasibinomial(link = 'logit')
  ),
  "Chose trans in private\noutcome round (=1)\n(pairs with trans only)" = glm(
    as.formula(
      str_glue("r2_choose_comparator ~ discussion_full + stratum_id + video_type_placebo + video_type_treatment + delivery_incentive_exp + {model_list[[3]]$lasso_controls %>% fml_sum()} + item_diff + r2_reliability_diff * r2_reliability_shown ")
    ),
    data = r2_choices_num_with_lasso %>% filter(discuss_type %in% c("control", "discussion_full")) %>% filter(!is.na(r2_choose_trans)),
    family = quasibinomial(link = 'logit')
  )
) %>%
  map(marginaleffects::avg_slopes, vcov = ~group_id)


logits_n_participants <- list("Num. participants" = logits %>%
  map(~ attr(.x, "model")$data) %>%
  map_int(~ n_distinct(.x$ind_id)) %>%
  as.character()
)

logits_n_groups <- list("Num. participants" = logits %>%
  map(~ attr(.x, "model")$data) %>%
  map_int(~ n_distinct(.x$group_id)) %>%
  as.character()
)


tex_export(
  logits,
  file = "outputs/tables/logits.tex",
  coef_rename = coef_label,
  gof_map = fe_label,
  coef_reorder = c("pair_includes_trans_discussion_full"),
  coef_omit = vars_to_regex("group_control|stratum_id|Intercept|video_type|delivery|pair_includes_trans_alt", control_vars),
  add_rows = combine_add_rows(
    list_to_add_rows(logits_n_participants),
    list_to_add_rows(logits_n_groups),
    list_to_add_rows(control_rows)
  )
)


# Long-run outcomes ---------------------------------------------------------------

# Does sensitivity to items change across treatments? ---------------------------------------------------------------

models_item_sensitivity <- list(
  # 1. item_diff weaker for trans
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff * pair_includes_trans + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video"
  ),

  #  2. item diff weaker for trans - doesn't change across treatments
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff * pair_includes_trans * discussion_full + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video"
  ),


  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff * discussion_full + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff_value_100 * discussion_full + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video"
  ),

  "Chose worker in\ntreatment round (=1)" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff * discussion_full + quality_diff * quality_included,
    data = r1_choices %>% filter(discuss_type %in% c("control", "discussion_full")),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Chose worker in\ntreatment round (=1)" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt * (video_type_placebo + video_type_treatment + discussion_full + factor(stratum_id) + delivery_incentive_exp) + item_diff_value_100 * discussion_full + quality_diff * quality_included,
    data = r1_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )
)

tex_export(
  models_item_sensitivity,
  file = "outputs/tables/item_sensitivity.tex",
  coef_rename = coef_label,
  coef_reorder = c(
    "pair_includes_trans_alt:discussion_full",
    "discussion_full",
    "item_diff",
    "discussion_full:item_diff"
  ),
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("pair_includes_trans_alt$", "stratum_id", "Intercept", "video", str_subset(control_vars, "item_diff", negate = TRUE), "delivery", "quality_diff")),
  add_rows = combine_add_rows(
    list_to_add_rows(list("Controls" = rep("X", length(models_item_sensitivity)))),
    list_to_add_rows(list("Controls interacted with worker is trans" = rep("X", length(models_item_sensitivity))))
  ),
  stat_vec = "long"
)

df


# Statistical discrimination ---------------------------------------------------------------


models_statistical_discrim <- list(
  "Chose worker in private\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ item_diff + pair_includes_trans * discussion_full + (video_type_placebo + video_type_treatment + factor(stratum_id) + delivery_incentive_exp) + (r2_reliability_diff * r2_reliability_shown ) * pair_includes_trans,
    data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    )
  ),
  "Chose worker in private\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ item_diff + pair_includes_trans * discussion_full + (video_type_placebo + video_type_treatment + factor(stratum_id) + delivery_incentive_exp) + (r2_reliability_diff * r2_reliability_shown ) * pair_includes_trans * discussion_full,
    data = r2_choices_with_item_val %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    coef_omit = "stratum_id|Intercept|video",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    )
  )
)

models_statistical_discrim[[2]] %>% tidy_90 %>% print_all


tex_export(
  models_statistical_discrim,
  file = "outputs/tables/statistical_discrim.tex",
  coef_rename = coef_label,
  coef_reorder = "pair_includes_trans:discussion_full",
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("group_control", "stratum_id", "Intercept", "video", str_subset(control_vars, "quality|reliability", negate = TRUE), "delivery_incentive")),
  controls_row = TRUE
)



models_statistical_discrim[[2]] %>% get_p_val("pair_includes_trans:discussion_full:r2_reliability_shown") %>%
  write_stat("outputs/stats/p_val_effect_on_statistical_discrim.tex")

models_statistical_discrim[[1]] %>% get_coeff("pair_includes_trans:r2_reliability_shown") %>%
  times_100 %>%
  write_stat("outputs/stats/statistical_discrim.tex", digits = 1)


# Divide by non-trans ---------------------------------------------------------------


# BY GENDER
r2_summ_mf <- r2_choices %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  group_by(discuss_type, pair_includes_female, pair_includes_trans) %>%
  summarise(mean_cl_boot(r2_choose_comparator),
            n_obs = n(),
            n_ind = n_distinct(ind_id),
            n_groups = n_distinct(group_id)
  ) %>%
  mutate(
    pair_type = case_when(
      pair_includes_trans ~ "Trans",
      pair_includes_female ~ "Female",
      TRUE ~ "Male"
    ),
    pair_type_label = pair_type %>%
      paste0(., "\n(Obs=", n_obs, ")")
  ) %>%
  mutate(
    discuss_label = ifelse(discuss_type=="discussion_full", "3-person discussion", "No discussion (private)") %>%
      paste0(., "\n(N=", n_ind, ")") %>%
      fct_relevel(c(first(sort(unique(.)))))
  )

r2_summ_mf %>%
  ggplot(aes(x = pair_type_label, y = y, ymin = ymin, ymax = ymax, fill = pair_type)) +
  geom_col(show.legend = FALSE, position = position_dodge(0.8)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.8)) +
  labs(y = "Probability of selecting alternative worker\nin outcome round", fill = element_blank(), x = "Gender of alternative worker\n(always compared to male benchmark)") +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
  geom_text(aes(y = 0.03, label = str_glue("{round(y, 3)*100}%")), position = position_dodge(0.8), size = 3) +
  facet_grid(~ discuss_label, scales = "free_x") +
  theme_minimal() +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(),
    axis.ticks = element_line(),
    panel.spacing = unit(0.7, "cm")
  )


ggsave("outputs/figs/r2_summ_mf.pdf", width = 6, height = 4)


r2_choices %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  filter(pair_includes_female) %>% pull(r2_choose_comparator) %>% mean_na() %>%
  write_percentage("outputs/stats/p_choose_female.tex")

r2_choices %>%
  filter(discuss_type %in% c("control")) %>%
  filter(pair_includes_trans) %>% pull(r2_choose_comparator) %>% mean_na() %>%
  write_percentage("outputs/stats/p_choose_trans_control.tex")

r2_choices %>%
  filter(discuss_type %in% c("discussion_full")) %>%
  filter(pair_includes_trans) %>% pull(r2_choose_comparator) %>% mean_na() %>%
  write_percentage("outputs/stats/p_choose_trans_treat.tex")

# Heterogeneity by gender PARTICIPANT
r2_choices %>%
  filter(discuss_type %in% c("discussion_full")) %>%
  filter(!pair_includes_trans) %>% pull(r2_choose_comparator) %>% mean_na() %>%
  write_percentage("outputs/stats/p_choose_non_trans_treat.tex")


# P-values - comparing trans to males
feols(
  r2_choose_comparator ~ pair_includes_trans,
  data = r2_choices %>% filter(discuss_type %in% c("control"), !pair_includes_female),
  cluster = "group_id"
) %>%
  get_p_val("pair_includes_transTRUE") %>%
  write_stat("outputs/stats/p_val_trans_vs_men_control.tex", digits = 2, p_value = TRUE)

feols(
  r2_choose_comparator ~ pair_includes_trans,
  data = r2_choices %>% filter(discuss_type %in% c("discussion_full"), !pair_includes_female),
  cluster = "group_id"
) %>%
  get_p_val("pair_includes_transTRUE") %>%
  write_stat("outputs/stats/p_val_trans_vs_men_treat.tex", digits = 2, p_value = TRUE)


r2_choices %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  group_by(discuss_type, female, pair_includes_female, pair_includes_trans) %>%
  summarise(mean_cl_cluster(r2_choose_comparator, cluster = group_id),
            n_obs = n(),
            n_ind = n_distinct(ind_id),
            n_groups = n_distinct(group_id)
  ) %>%
  rename(y = estimate,
         ymin = conf.low, ymax = conf.high) %>%
  mutate(
    pair_type = case_when(
      pair_includes_trans ~ "Trans",
      pair_includes_female ~ "Female",
      TRUE ~ "Male"
    ),
    pair_type_label = pair_type
  ) %>%
  mutate(
    female_label = ifelse(female == 1, "Female participant", "Male participant")
  ) %>%
  mutate(
    discuss_label = ifelse(discuss_type=="discussion_full", "3-person discussion", "No discussion (private)") %>%
      fct_relevel(c(first(sort(unique(.)))))
  ) %>%
  ggplot(aes(x = pair_type_label, y = y, ymin = ymin, ymax = ymax, fill = pair_type)) +
  geom_col(show.legend = FALSE, position = position_dodge(0.8)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.8)) +
  labs(y = "Probability of selecting alternative worker\nin outcome round", fill = element_blank(), x = "Gender of alternative worker\n(always compared to male benchmark)") +
  scale_y_continuous(labels = scales::percent, expand = c(0.01, 0), limits = c(0, 0.8)) +
  geom_text(aes(y = 0.03, label = str_glue("{round(y, 3)*100}%")), position = position_dodge(0.8), size = 3) +
  facet_grid(female_label ~ discuss_label, scales = "free_x") +
  theme_bw() +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major = element_blank()
  )

ggsave("outputs/figs/r2_by_participant_gender.pdf", width = 6, height = 6)


feols_custom(
  r2_choose_comparator ~ discussion_full,
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "phase"),
  cluster = "group_id",
  data = r2_choices %>% filter(discuss_type %in% c("discussion_full", "control"), pair_includes_female, female == 0),
  ri = FALSE,
  lasso = TRUE,
  lasso_options = list(
    potential_controls = control_vars,
    t = "discussion_full",
    interact = NULL,
    group_control = TRUE
  )
)

r2_choices$female


# Correlation between leadership scale and pro-trans in control group ---------------------------------------------------------------

ls %>% dups_report(ls_ind_id)

ls_scores_ind %>% dups_report(ls_ind_id)

feols_custom(
  r2_choose_trans ~ ls_score_fact_z,
  data = r2_with_ls %>% filter(discuss_type == "control", phase == "phase_2"),
  fixef = c("stratum_id", "delivery_incentive_exp", "video_type"),
  cluster = "group_id",
  ri = FALSE,
  n_sims = ri_sims,
  stratum = "stratum_id",
  var_types = var_types_spec,
  coef_omit = "stratum_id|Intercept|video"
)

# Discussion efffect - for those who refused audio consent ---------------------------------------------------------------

r2_choices_num %>%
  count_prop(discuss_type, audio_refused)


models_audio_consent <- list(
  "Chose trans in\noutcome round (=1)\n(pairs with trans only)" = feols_custom(
    r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_full + discussion_full:audio_refused + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full"), !is.na(r2_choose_trans)),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = NULL,
      group_control = TRUE
    )
  )
)


prop_audio_refused <- models_audio_consent[[1]]$data %>%
  filter(discussion_full == 1) %>%
  summarise(mean_na(audio_refused))

tex_export(
  models_audio_consent,
  file = "outputs/tables/audio_consent_robustness.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id", "Intercept", "video", "reliability", str_subset(control_vars, "r2_reliability", negate = TRUE), "delivery")),
  add_rows = tibble(
    l = "Mean: Audio recording refused $\\vert$ 3-person discussion",
    m1 = round(prop_audio_refused, 2)
  ),
  controls_row = TRUE
)


models_audio_consent[[1]] %>% get_p_val("discussion_full:audio_refusedTRUE") %>%
  write_stat("outputs/stats/p_val_audio_refused.tex")

r2_choices_num %>% bar_chart(
  x = discuss_type,
  y = r2_choose_trans,
  fill = audio_refused
)


# Correlation between own choices and group predic misperception ---------------------------------------------------------------

# Get individual-level dataset of r2 choices (number of times selected trans)
r2_ind <- r2_choices_num %>%
  group_by(across(all_of(id_vars))) %>%
  summarise(n_r2_choose_trans = sum_na(r2_choose_trans)) %>%
  ungroup


r2_ind %>% hist_basic(n_r2_choose_trans)

# Merge r2 choices onto the group_predic data set
group_predic_w_own_r2_choices <- group_predic_w_actual_choices %>%
  select(ind_id, round, group_predic_member_id, group_predic_includes_trans,
         group_predic_choose_trans, r2_choose_trans_other = r2_choose_trans,
         misper) %>%
  tidylog::left_join(r2_ind, by = "ind_id")


# Misperception correlation
group_predic_w_own_r2_choices %>%
  filter(discuss_type == "control") %>%
  bar_chart(x = n_r2_choose_trans, y = misper)

# Correlation between own choices and predicted others choices
group_predic_w_own_r2_choices %>%
  filter(discuss_type == "control") %>%
  bar_chart(x = n_r2_choose_trans, y = group_predic_choose_trans)

group_predic_long %>%
  tidylog::left_join(r2_ind, by = "ind_id", suffix = c("", "_r2")) %>%
  filter(discuss_type == "control" | discuss_type == "discussion_full") %>%
  bar_chart(x = n_r2_choose_trans, y = value, fill = name,
            facet = discuss_type)

# If social image is a function of (predic - own choices), then:
# - discriminators (n=0) should discriminate less
# - middle people (n=1) should discriminate slightly more
# - pro-trans people (n=2) should discriminate more

# Implies that social image would lead to compression of the distribution

# Does social image lead to compression in R1? ---------------------------------------------------------------

r1_choices_ind

icc_r1 <- r1_choices_ind %>%
  left_join(df %>% select(ind_id, phase, discuss_type), by = "ind_id") %>%
  group_by(discuss_type, phase) %>%
  nest() %>%
  rowwise() %>%
  mutate(
    icc = list(as_tibble(ICCest(x = group_id, y = p_self_selected_trans, data = data)))
  ) %>%
  unnest(icc)

# Higher ICC in R1 in public choices than in control
# People manage to coordinate and match each others' choices in social image...
# Everyone always seems the same choices...

r2_choices

# And R2?
icc_r2 <- r2_choices %>%
  group_by(ind_id, group_id, phase, discuss_type) %>%
  summarise(n_r2_choose_trans = sum_na(r2_choose_trans)) %>%
  group_by(discuss_type, phase) %>%
  nest() %>%
  rowwise() %>%
  mutate(
    icc = list(as_tibble(ICCest(x = group_id, y = n_r2_choose_trans, data = data)))
  ) %>%
  unnest(icc)

# Both phase 1 and 2
icc_r2_pooled <- r2_choices %>%
  group_by(ind_id, group_id, phase, discuss_type) %>%
  summarise(n_r2_choose_trans = sum_na(r2_choose_trans)) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  group_by(discuss_type) %>%
  nest() %>% rowwise() %>%
  mutate(
    icc = list(as_tibble(ICCest(x = group_id, y = n_r2_choose_trans, data = data)))
  ) %>%
  unnest(icc)


bind_rows(
  r1 = icc_r1,
  r2 = icc_r2,
  .id = "round"
) %>%
  ggplot(aes(x = discuss_type, y = ICC, ymin = LowerCI, ymax = UpperCI, colour = phase)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.5)) +
  geom_point(position = position_dodge(0.5)) +
  facet_wrap(~round)


# RI p-value of control vs choosing only
diff_icc <- function(df) {
  if (!("y" %in% names(df))) {
    stop("y must be specified")
  }

  suppressWarnings(
    icc_choosing_only <- df %>%
      filter(treat == 1) %>%
      ICCbare(x = group_id, y = y, data = .)
  )


  suppressWarnings(
    icc_control <- df %>%
      filter(treat == 0) %>%
      ICCbare(x = group_id, y = y, data = .)
  )


  diff <- icc_choosing_only - icc_control
  return(diff)
}

r1_for_ri <- r1_choices_ind %>%
  left_join(df %>% select(ind_id, phase, discuss_type, stratum_id), by = "ind_id") %>%
  filter(phase == "phase_2", (discuss_type == "choosing_only" | discuss_type == "control")) %>%
  mutate(y = as.numeric(p_self_selected_trans)) %>%
  mutate(choosing_only = as.numeric(discuss_type == "choosing_only"),
         treat = choosing_only) %>% 
  ungroup

r2_for_ri <- r2_choices %>%
  group_by(ind_id, group_id, phase, discuss_type, stratum_id) %>%
  summarise(n_r2_choose_trans = sum_na(r2_choose_trans)) %>%
  filter(phase == "phase_2", (discuss_type == "choosing_only" | discuss_type == "control")) %>%
  mutate(y = as.numeric(n_r2_choose_trans)) %>%
  mutate(choosing_only = as.numeric(discuss_type == "choosing_only"),
         treat = choosing_only) %>% 
  ungroup

r2_for_ri_pooled <- r2_choices %>%
  group_by(ind_id, group_id, phase, discuss_type, stratum_id) %>%
  summarise(n_r2_choose_trans = sum_na(r2_choose_trans)) %>%
  filter(phase == "phase_2", (discuss_type == "discussion_full" | discuss_type == "control")) %>%
  mutate(y = as.numeric(n_r2_choose_trans)) %>%
  mutate(discussion_full = as.numeric(discuss_type == "discussion_full"),
         treat = discussion_full) %>% 
  ungroup

ri_sims <- 500
set.seed(12345)

ri_icc <- list(
  r1 = conduct_ri_custom(
    test_function = diff_icc,
    declarations = list(
      get_ra_declaration(
        df = r1_for_ri,
        treat = "treat",
        stratum = "stratum_id",
        cluster = "group_id"
      )
    ),
    p = "two-tailed",
    data = r1_for_ri,
    sims = ri_sims,
    progress_bar = TRUE,
    assignments = list("treat")
  ),
  r2 = conduct_ri_custom(
    test_function = diff_icc,
    declarations = list(
      get_ra_declaration(
        df = r2_for_ri,
        treat = "treat",
        stratum = "stratum_id",
        cluster = "group_id"
      )
    ),
    p = "two-tailed",
    data = r2_for_ri,
    sims = ri_sims,
    progress_bar = TRUE,
    assignments = list("treat")
  ),
  r2_pooled = conduct_ri_custom(
    test_function = diff_icc,
    declarations = list(
      get_ra_declaration(
        df = r2_for_ri_pooled,
        treat = "treat",
        stratum = "stratum_id",
        cluster = "group_id"
      )
    ),
    p = "two-tailed",
    data = r2_for_ri_pooled,
    sims = ri_sims,
    progress_bar = TRUE,
    assignments = list("treat")
  )
)

summary.ri <- ri2:::summary.ri


# R1 - phase 2 - choosing only vs control
# R2 - phase 2 - choosing only vs control
# R2 - Phase 1 + 2 - control vs discussion full


icc_table <- bind_rows(
  r1 = icc_r1 %>% filter(phase == "phase_2", discuss_type %in% c("control", "choosing_only")),
  r2 = icc_r2 %>% filter(phase == "phase_2", discuss_type %in% c("control", "choosing_only")),
  r2 = icc_r2_pooled %>% mutate(phase = "Phases 1 + 2"),
  .id = "type"
) %>% 
  mutate(
    type = recode(
      type,
      r1 = "Treatment",
      r2 = "Outcome"
    ),
    phase = recode(
      phase,
      "phase_2" = "Phase 2 only",
      "Phases 1 + 2" = "Phases 1 + 2"
    ),
    discuss_type = recode(
      discuss_type,
      "control" = "(a) No discussion (private)",
      "choosing_only" = "(b) No discussion (public)",
      "discussion_full" = "(b) 3-person discussion"
    )
  ) %>%
  mutate(ci = paste0("[", format(round(LowerCI, 2), nsmall = 2), ", ", format(round(UpperCI, 2), nsmall = 2), "]")) %>%
  select(Treatment = discuss_type, Sample = phase,  Round = type, ICC = ICC, `95% CI` = ci, `N groups` = N)

icc_table$Round[c(2, 4, 6)] <- ""
icc_table$Sample[c(2, 4, 6)] <- ""
icc_table$`p-val: (a)=(b)` <- rep(NA, 6)
icc_table$`p-val: (a)=(b)`[c(1, 3, 5)] <- c(
  ri_icc %>%
    map(ri2:::summary.ri) %>%
    map_dbl("two_tailed_p_value")
)

icc_table %>%
  xtable(digits = 2) %>%
  print(file = "outputs/tables/icc.tex", floating = FALSE, include.rownames = FALSE,
        hline.after = c(-1, 0, 2, 4, 6))


icc_table$`p-val: (a)=(b)`[[1]] %>%
  write_stat("outputs/stats/ri_p_icc_private_public.tex")


# For magnitudes intuition, plot the correlation between R1 choices with others for control and public
r1_choices %>%
  group_by(group_id, ind_id, discuss_type) %>%
  summarise(n_trans_ind = sum_na(r1_choose_trans),
            p_trans_ind = mean_na(r1_choose_trans)) %>%
  group_by(group_id) %>%
  mutate(n_trans_group = sum_na(n_trans_ind)) %>%
  ungroup %>%
  mutate(n_trans_other = n_trans_group - n_trans_ind,
         p_trans_other = n_trans_other / 4) %>%
  mutate(p_trans_other = round(p_trans_other, 2)) %>%

  filter(discuss_type == "control" | discuss_type == "choosing_only") %>%
  group_by(discuss_type, p_trans_other) %>%
  summarise(mean_ci(p_trans_ind)) %>%

  ggplot(aes(x = p_trans_other, y = y, ymin = ymin, ymax = ymax, colour = discuss_type), alpha = 0.7) +
  geom_point(position = position_dodge(0.05)) + geom_errorbar(position = position_dodge(0.05), width = 0.02) +
  labs(x = "P(others selected trans in your group)", y = "P(selecting trans yourself)") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  coord_cartesian(xlim = c(0, 1), ylim = c(0, 1))



# x = others' choices
# y = my choices


# ICC by relation scores ----------------------------------------------------------------------------------------------

icc_r1_by_relation <- r1_choices_ind %>%
  left_join(df %>% select(ind_id, discuss_type, stratum_id), by = "ind_id") %>%
  left_join(relations_scores_group, by = "group_id") %>%
  mutate(relation_score_fact_z_group = relation_score_fact_z_group > median_na(relation_score_fact_z_group)) %>%
  group_by(discuss_type, relation_score_fact_z_group) %>%
  nest() %>% 
  rowwise() %>%
  mutate(
    icc = list(as_tibble(ICCest(x = group_id, y = p_self_selected_trans, data = data)))
  ) %>%
  unnest(icc) %>%
  arrange(discuss_type, relation_score_fact_z_group)

# Create plot of ICC by relation score
icc_r1_by_relation %>%
  mutate(
    relation_score_fact_z_group = if_else(relation_score_fact_z_group, "Above median", "Below median"),
    relation_score_fact_z_group = factor(relation_score_fact_z_group, levels = c("Below median", "Above median"))
  ) %>%
  ggplot(aes(x = discuss_type, y = ICC, ymin = LowerCI, ymax = UpperCI, colour = relation_score_fact_z_group)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.5)) +
  geom_point(position = position_dodge(0.5)) +
  labs(
    x = "Discussion type",
    y = "Intraclass correlation coefficient",
    colour = "Group relation score"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")


  # FOCUS ON CONTROL GROUP
  r1_choices_ind %>%
  left_join(df %>% select(ind_id, discuss_type, stratum_id), by = "ind_id") %>%
  left_join(relations_scores_group, by = "group_id") %>% 
  filter(discuss_type == "control" | discuss_type == "choosing_only") %>%
  feols_custom(p_self_selected_trans ~ p_other_selected_trans_control * relation_score_fact_z_group, data = ., cluster = "group_id", fixef = "stratum_id")


# ------------------------------------------------------------------------------
# New Randomization Inference: Compare ICC in the control group
# for groups with close relations vs. groups with not close relations.
# We restrict to discuss_type == "control", use the outcome round (r2 choices)
# and create a binary treat variable using the relation scores.
# ------------------------------------------------------------------------------

# Create a new data frame for RI analysis by filtering only the control group 
# and joining with the group relation scores.
r2_for_ri_relation <- r2_for_ri %>%
  filter(discuss_type == "control") %>%
  left_join(relations_scores_group, by = "group_id") %>%
  # Create binary indicator: 1 if group relation score is above the overall median,
  # otherwise 0.
  ungroup %>% 
  mutate(y = as.numeric(n_r2_choose_trans)) %>%
  mutate(treat = as.numeric(relation_score_fact_z_group > median(relation_score_fact_z_group, na.rm = TRUE))) %>% 
  ungroup

  r2_for_ri_relation %>% count_prop(treat)

ri_sims <- 500
set.seed(12345)

# Run randomization inference using diff_icc.
ri_icc_relation <- conduct_ri_custom(
  test_function = diff_icc,
  declarations = list(
    get_ra_declaration(
      df = r2_for_ri_relation,
      treat = "treat",
      stratum = "stratum_id",
      cluster = "group_id"
    )
  ),
  p = "two-tailed",
  data = r2_for_ri_relation,
  sims = ri_sims,
  progress_bar = TRUE,
  assignments = list("treat")
)

# Print a summary of the RI results.
print(summary.ri(ri_icc_relation))

# Optionally, write the two-tailed p-value of the RI result to file.
ri_icc_pval <- summary.ri(ri_icc_relation)$two_tailed_p_value
write_stat(ri_icc_pval, "outputs/stats/ri_p_icc_relation_control.tex", digits = 2, p_value = TRUE)


# Inferring the relative persuasiveness of discussion vs observing choices ---------------------------------------------------------------

# I want to compare choosing only observers to listeners...
# and CONDITION on the proportion of pro-trans choices they observed in first round

r2_persuasion <- r2_with_announce %>%
  filter(
    observer == 1 | is_listener == 1 | discuss_type == "control"
  ) %>%
  mutate(is_listener = ifelse(is.na(is_listener), 0, is_listener)) %>%
  count_prop(discuss_type, observer, is_listener) %>%
  mutate(saw_choices = observer == 1 | is_listener == 1) %>%
  mutate(heard_reasons = is_listener == 1)


r2_persuasion %>%
  filter(discuss_type != "control") %>%
  mutate(is_listener = as.logical(is_listener)) %>%
  bar_chart(x = p_other_selected_trans, y = r2_choose_trans, fill = is_listener, return_data = TRUE) %>%
  ggplot(aes(x = p_other_selected_trans, y = y, ymin= ymin, ymax = ymax, colour = is_listener)) +
  geom_point() +
  geom_errorbar()



models_persuasion <- list(
  feols_custom(
    r2_choose_trans ~ saw_choices + heard_reasons + p_other_selected_trans,
    data = r2_persuasion,
    fixef = c("stratum_id", "delivery_incentive_exp"),
    cluster = "group_id"
  ),
  feols_custom(
    r2_choose_trans ~ saw_choices + heard_reasons + p_other_selected_trans + saw_choices*p_other_selected_trans + heard_reasons * p_other_selected_trans,
    data = r2_persuasion,
    fixef = c("stratum_id", "delivery_incentive_exp"),
    cluster = "group_id"
  )
)

# Is the gradient with respect to p_other_selected_trans greater for listeners than for observers??
# At the moment no


# SDB robustness ---------------------------------------------------------------

sdb_robustness <- list(
  "SDB:\nacquiescence bias\ncorrection" = feols_custom(
    r2_choose_trans ~ discussion_full * high_sdb_acqui + video_type  + factor(stratum_id),
    data = r2_with_sdb %>% filter(!is.na(r2_choose_trans)) %>% filter(discuss_type %in% c("discussion_full", "control"), phase == "phase_1"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = NULL,
      group_control = TRUE
    )
  ),
  "SDB:\nfactor\nanalysis" = feols_custom(
    r2_choose_trans ~ discussion_full * high_sdb_fact + video_type  + factor(stratum_id),
    data = r2_with_sdb %>% filter(!is.na(r2_choose_trans)) %>% filter(discuss_type %in% c("discussion_full", "control"), phase == "phase_1"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = NULL,
      group_control = TRUE
    )
  ),
  "SDB:\ninverse covariance\nweighted" = feols_custom(
    r2_choose_trans ~ discussion_full * high_sdb_inv_cov + video_type  + factor(stratum_id),
    data = r2_with_sdb %>% filter(!is.na(r2_choose_trans)) %>% filter(discuss_type %in% c("discussion_full", "control"), phase == "phase_1"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = NULL,
      group_control = TRUE
    )
  )
)

tex_export(
  sdb_robustness,
  file = "outputs/tables/sdb_robustness.tex",
  coef_omit = vars_to_regex("stratum_id|Intercept|video_type|items", control_vars),
  coef_rename = coef_label,
  controls_row = TRUE,
  gof_map = fe_label_no_fe,
  additional_header = vec_to_custom_header(
    c(" ", rep("Dep var: Chose trans in private outcome round (=1)", 3))
  )
)


# High-stakes condition ---------------------------------------------------------------


ggred <- gg_color_hue(2)[[1]]
ggblue <- gg_color_hue(2)[[2]]

r2_summ_delivery_incentive <- r2_choices %>%
  filter(!is.na(delivery_incentive)) %>%
  group_by(group_label, discuss, discuss_type_label, discuss_type, delivery_incentive, pair_includes_trans) %>%
  summarise(mean_cl_boot(r2_choose_comparator),
            n_obs = n(),
            n_ind = n_distinct(ind_id),
            n_groups = n_distinct(group_id)
  ) %>%
  mutate(
    discuss_label = ifelse(discuss, "Discussion", "Control"),
    pair_includes_trans_label = ifelse(pair_includes_trans, "Worker\nis trans", "Worker\nis non-trans")
  ) %>%
  mutate(delivery_incentive = fct_relevel(delivery_incentive, "multi", "single")) %>%

  mutate(
    trans_and_delivery = paste0(pair_includes_trans, "_", delivery_incentive) %>%
      factor(),
    trans_and_delivery = fct_relevel(trans_and_delivery, "FALSE_single", "FALSE_multi", "TRUE_single", "TRUE_multi")
  )


r2_summ_delivery_incentive %>%

  ggplot(aes(x = pair_includes_trans_label, y = y, ymin = ymin, ymax = ymax, fill = trans_and_delivery)) +
  geom_col(aes(colour = trans_and_delivery), width = 0.8, position = position_dodge(0.8)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.8), show.legend = FALSE) +
  labs(y = "Probability of selecting the worker", fill = element_blank(), x = element_blank()) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0), limits = c(0, 0.76)) +
  geom_text(aes(y = 0.03, label = str_glue("{round(y, 3)*100}%")), show.legend = FALSE, position = position_dodge(0.8), size = 3) +
  facet_wrap(~ discuss_type_label, scales = "free_x") +
  theme_minimal() +
  scale_fill_manual(
    name = "",
    values = c(ggred, "white", ggblue, "white"),
    labels = c("Worker is non-trans, 1 delivery", "Worker is non-trans, 3 deliveries",
               "Worker is T, 1 delivery", "Worker is T, 3 deliveries")
  ) +
  scale_color_manual(
    name = "",
    values = c(ggred, ggred, ggblue, ggblue),
    labels = c("Worker is non-trans, 1 delivery", "Worker is non-trans, 3 deliveries",
               "Worker is T, 1 delivery", "Worker is T, 3 deliveries")
  ) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    panel.grid.major = element_blank(),
    axis.line = element_line(),
    axis.ticks = element_line(),
    panel.spacing = unit(0.7, "cm")
  )


df %>% count_prop(delivery_incentive_exp)

ggsave("outputs/figs/delivery_incentive.pdf", width = 9, height = 6, scale = 1)


# Is treatment effect stronger in multi delivery condition?
feols_custom(
  r2_choose_comparator ~ group_label + delivery_incentive + group_label * delivery_incentive,
  data = r2_choices %>% filter(pair_includes_trans, !is.na(delivery_incentive)),
  fixef = c("stratum_id", "video_type"),
  cluster = "group_id"
)


feols_custom(
  r2_choose_comparator ~ group_label,
  data = r2_choices %>% filter(pair_includes_trans, delivery_incentive == "multi"),
  fixef = c("stratum_id", "video_type"),
  cluster = "group_id"
)

r2_choices %>% count_prop(delivery_incentive_exp)

models_high_stakes <- list(
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * (discussion_full) * delivery_incentive_multi,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full"))%>% filter(!is.na(delivery_incentive)),
    fixef = NULL,
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id|delivery_incentive_exp",
  ),
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt*(video_type_placebo + video_type_treatment + discussion_full * delivery_incentive_multi + factor(stratum_id) + delivery_incentive_multi) + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase)%>% filter(!is.na(delivery_incentive)),
    fixef = c("stratum_id", "video_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id|delivery_incentive_exp"
  ),
  "Chose trans in\noutcome round (=1)\n(pairs with trans only)" = feols_custom(
    r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_full* delivery_incentive_multi + factor(stratum_id) + delivery_incentive_multi + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% filter(!is.na(r2_choose_trans))%>% filter(!is.na(delivery_incentive)),
    fixef = c("stratum_id", "video_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id"
  )
)

tex_export(
  models_high_stakes,
  file = "outputs/tables/high_stakes.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id", "Intercept", "video", "group_control", "pair_includes_trans_alt$", control_vars)),
  add_rows = list_to_add_rows(control_rows)
)

models_high_stakes[[1]] %>% get_p_val("pair_includes_trans:discussion_full:delivery_incentive_multi") %>%
  write_stat("outputs/stats/high_stakes_p_val.tex", 2, p_value = TRUE)

n_high_stakes_sample <- df %>% filter(!is.na(delivery_incentive)) %>%
  nrow()

n_high_stakes_sample %>% write_stat("outputs/stats/n_high_stakes_sample.tex", 0, p_value = FALSE)

n_high_stakes <- df %>% count_prop(delivery_incentive) %>%
  filter(delivery_incentive == "multi") %>%
  nrow()

n_high_stakes %>% write_stat("outputs/stats/n_high_stakes.tex", 0, p_value = FALSE)

n_low_stakes <- (n_high_stakes_sample - n_high_stakes)

n_low_stakes %>% write_stat("outputs/stats/n_low_stakes.tex", 0, p_value = FALSE)


feols_custom(
  r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt*(video_type_placebo + video_type_treatment + discussion_full * delivery_incentive_multi + factor(stratum_id) + delivery_incentive_multi) + item_diff + r2_reliability_diff * r2_reliability_shown ,
  data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase)%>% filter(!is.na(delivery_incentive)),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id",
  ri = FALSE,
  n_sims = ri_sims,
  stratum = "stratum_id",
  lasso = TRUE,
  lasso_options = list(
    potential_controls = control_vars,
    t = "discussion_full",
    interact = "pair_includes_trans",
    group_control = TRUE
  ),
  var_types = var_types_spec,
  coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
  coef_omit = "group_control|stratum_id|delivery_incentive_exp"
) %>% tidy_90 


df


# High stakes - video ---------------------------------------------------------------

models_high_stakes_video <- list(
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * (video_type_placebo + video_type_treatment) * delivery_incentive_multi,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full"))%>% filter(!is.na(delivery_incentive)),
    fixef = NULL,
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id|delivery_incentive_exp",
  ),
  "Chose worker in\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt*(discussion_full + (video_type_placebo + video_type_treatment) * delivery_incentive_multi + factor(stratum_id) + delivery_incentive_multi) + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% count_prop(phase)%>% filter(!is.na(delivery_incentive)),
    fixef = c("stratum_id", "discuss_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id|delivery_incentive_exp"
  ),
  "Chose trans in\noutcome round (=1)\n(pairs with trans only)" = feols_custom(
    r2_choose_trans ~ discussion_full + (video_type_placebo + video_type_treatment) * delivery_incentive_multi + factor(stratum_id) + delivery_incentive_multi + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% filter(!is.na(r2_choose_trans))%>% filter(!is.na(delivery_incentive)),
    fixef = c("stratum_id", "discuss_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = "discussion_full",
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    var_types = var_types_spec,
    coef_keep = "^pair_includes_trans$|(^pair_includes_trans\\:discussion_full$)|^discussion_full$|item_diff|r2_reliability_diff|r2_reliability_shown",
    coef_omit = "group_control|stratum_id"
  )
)

tex_export(
  models_high_stakes_video,
  file = "outputs/tables/high_stakes_video.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id", "Intercept", "discussion", "group_control", "pair_includes_trans_alt$", control_vars)),
  add_rows = list_to_add_rows(control_rows)
)

r2_choices_num %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  filter(!is.na(r2_choose_trans))%>% filter(!is.na(delivery_incentive)) %>%
  bar_chart(
    x = video_type,
    y = r2_choose_trans,
    fill = delivery_incentive
  )


# Context  ---------------------------------------------------------------


# Familiarity of deliveries ---------------------------------------------------------------

df %>%
  mutate(across(c(i1, i2, i3), \(x) ifelse(x == -98, 0, x))) %>%
  mutate(market_research = b11 | b12,
         hired_someone = b8,
         app_order = i1 | i2 | i3) %>%
  pivot_longer(c(hired_someone, market_research, app_order, b11)) %>%
  bar_chart(x = name, y = value, fill = name, flip = FALSE, perc = TRUE, val_label = TRUE)

familiarity <- df %>%
  mutate(market_research = b11 | b12,
         hired_someone = b8,
         app_order = i1 | i2 | i3)


familiarity %>%
  pull(app_order) %>% mean_na() %>%
  write_percentage("outputs/stats/have_ordered_from_app.tex")


familiarity %>%
  pull(market_research) %>% mean_na() %>%
  write_percentage("outputs/stats/taken_part_market_research.tex")

# Comprehension checks ---------------------------------------------------------------

practice_check <- df %>%
  select(matches("practice_check")) %>%
  mutate(
    practice_check_2_correct = practice_check_2 == 4,
    practice_check_3_correct = practice_check_3 == 1,
    practice_check_7_correct = practice_check_7 == 1,
    practice_check_4_correct = practice_check_4 == 1
  )

practice_check %>%
  pivot_longer(matches("practice_check_\\d_correct")) %>%
  bar_chart(x = name, y = value, flip = TRUE, val_label = TRUE)




practice_check %>%
  pivot_longer(matches("practice_check_\\d_correct")) %>%
  pull(value) %>%
  mean_na() %>%
  write_percentage("outputs/stats/practice_check_correct.tex")

hiring_check <- df %>%
  select(matches("hiring_check"), matches("delivery_incentive")) %>%
  mutate(
    hiring_check_1_correct = hiring_check_1 == "1 2 3",
    hiring_check_5_correct = hiring_check_5 == 10,
    hiring_check_3_correct = ifelse(
      delivery_incentive == "multi" & !is.na(delivery_incentive),
      hiring_check_3 == 3,
      hiring_check_3 == 1
    ),
    hiring_check_5.1_correct = hiring_check_5.1 == 1,
    hiring_check_11_correct = hiring_check_11 == 1
  )


hiring_check %>%
  pivot_longer(matches("hiring_check_(.*)_correct")) %>%
  bar_chart(x = name, y = value, flip = TRUE, val_label = TRUE)

hiring_check %>%
  pivot_longer(matches("hiring_check_(.*)_correct")) %>%
  pull(value) %>%
  mean_na() %>%
  write_percentage("outputs/stats/hiring_check_correct.tex")


# Order effects ---------------------------------------------------------------


# P-value of order effects for treatment round

order_effect_models <- list(
  "alt_model" = feols(
    choose_comparator ~ round_2 + round_3 + round_4 + discussion_full * (round_2 + round_3 + round_4),
    order_effects_df %>% filter(discuss_type %in% c("control", "discussion_full"), as.numeric(pair_includes_trans)==1),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id"
  ))


p_value_order_effects <- f_test(order_effect_models$alt_model, "round_\\d\\:discussion_full")$p

p_value_order_effects %>%
  write_stat("outputs/tables/p_value_order_effects.tex", digits = 2)


p_lab <- paste0(round(p_value_order_effects, 2))
p_lab_start <- bquote("F-test ("*beta[r1]~"="*beta[r2]~"="*beta[r3]~"="*beta[r4]*"):") %>% as.expression()
p_lab_second <- paste0("p=", p_lab)

plot_order <- order_effects_df %>%
  select(choose_comparator, discuss_type, discuss_type_label, trans_round, round_type) %>%
  group_by(discuss_type, discuss_type_label, trans_round, round_type) %>%
  summarise(
    mean_cl_boot(choose_comparator),
    n_obs = n()
  ) %>%

  filter(discuss_type %in% c("control", "discussion_full")) %>%

  ggplot(aes(x = trans_round, y = y, ymin = ymin, ymax = ymax, fill = discuss_type_label)) +
  geom_col(width = 0.8, position = position_dodge(0.8)) +
  geom_errorbar(width = 0.2, position = position_dodge(0.8)) +
  facet_wrap(~ round_type, scales = "free_x") +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0), lim = c(0, 0.85)) +
  labs(x = "Order of delivery option pair that includes transgender", y = "Proportion choosing transgender",
       fill = element_blank()) +
  theme_bw() +
  scale_color_brewer(palette = "Set2", direction = -1) +
  scale_fill_brewer(palette = "Set2", direction = -1) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())


plot_order %>% annotate_figure(
  top = text_grob(p_lab_start, hjust = -0.8, vjust = 5, size = 10)
) %>%
  annotate_figure(
    top = text_grob(p_lab_second, hjust = -2.9, vjust = 10.8, size = 10)
  )

ggsave("outputs/figs/order_effects.pdf", width = 6, height = 4)


# Increased "polarisation" in discussion groups? ---------------------------------------------------------------

trans_per_group <- r2_choices %>%
  filter(pair_includes_trans) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  group_by(group_id, group_label, discuss_type, discuss_type_label) %>%
  summarise(n_choose_trans = sum(r2_choose_trans, na.rm = TRUE))


cdf_basic <- trans_per_group %>%
  group_by(discuss_type_label, n_choose_trans) %>%
  summarise(n = n()) %>%
  arrange(discuss_type_label, n_choose_trans) %>%
  group_by(discuss_type_label) %>%
  mutate(cdf = cumsum(n) / sum(n))

cdf_dashed <- cdf_basic %>%
  mutate(type = "cdf") %>%
  bind_rows(cdf_basic %>%
              mutate(type = "prior",
                     cdf = lag(cdf))) %>%
  drop_na() %>%
  ungroup %>%
  add_row(
    discuss_type_label = "No discussion (private)",
    n_choose_trans = 0,
    n = NA,
    cdf = NA,
    type = "prior"
  ) %>%
  arrange(discuss_type_label, n_choose_trans, desc(type))

cdf_dashed

cdf_dashed %>%
  group_by(discuss_type_label) %>%
  ggplot(aes(x = n_choose_trans, y = cdf, colour = discuss_type_label)) +
  geom_point(aes(colour = discuss_type_label, shape = type)) +
  scale_shape_manual(values = c(16, 21)) +
  geom_segment(
    aes(x = lag(n_choose_trans), y = lag(cdf), colour = discuss_type_label,
        xend = n_choose_trans, yend = cdf,
        lty = fct_rev(type))) +
  theme_classic() +
  labs(x = "Number of times trans is selected by group in outcome round (/6)",
       y = "Cumulative proportion of groups",
       colour = element_blank()) +
  scale_color_brewer(palette = "Set2", direction = 1) +
  scale_y_continuous(breaks = seq(0, 1, 0.2), expand = c(0, 0), lim = c(0, 1.1)) +
  guides(linetype = FALSE, shape = FALSE)


ggsave("outputs/figs/cdf_choose_trans_in_group.pdf", width = 6, height = 4)


# More or less sensitive to other characteristics? ---------------------------------------------------------------


sensitivity_group <- feols_custom(
  r2_choose_comparator ~ pair_includes_trans * (discussion_full + factor(stratum_id)) + r2_reliability_diff * r2_reliability_shown  + delivery_incentive_exp + factor(item_diff),
  data = r2_choices %>% filter(discuss_type == "discussion_full"),
  fixef = c("stratum_id", "comparator_order_in_pair"),
  cluster = "group_id",
)

f_test(sensitivity_group, "r2_reliability|item_diff")


sensitivity_ind <- feols_custom(
  r2_choose_comparator ~ pair_includes_trans * (group_label + factor(stratum_id)) + r2_reliability_diff * r2_reliability_shown  + delivery_incentive_exp + factor(item_diff),
  data = r2_choices %>% filter(discuss_type == "control"),
  fixef = c("stratum_id", "comparator_order_in_pair"),
  cluster = "group_id",
)

f_test(sensitivity_ind, "r2_reliability|item_diff")


# Get RI declaration for both discussion full and pair_includes_trans
ri_discussion_full <- get_ra_declaration(
  df = r2_choices %>% filter(discuss_type %in% c("control", "discussion_full")),
  treat = "discussion_full",
  stratum = "stratum_id",
  cluster = "group_id"
)

ri_pair_includes_trans <- get_ra_declaration(
  df = r2_choices %>% filter(discuss_type %in% c("control", "discussion_full")),
  treat = "pair_includes_trans",
  within = "ind_id"
)



# GET p-value for difference between f-stats
compare_f_stats <- function(df) {
  f1 <- df %>% filter(discussion_full == 0) %>%
    feols_custom(
      r2_choose_trans ~ factor(stratum_id) + r2_reliability_diff * r2_reliability_shown  + delivery_incentive_exp + factor(item_diff),
      data = .,
      fixef = c("stratum_id", "comparator_order_in_pair"),
      cluster = "group_id",
    ) %>%
    f_test("r2_reliability|item_diff", print = FALSE) %>% .$stat

  f2 <- df %>% filter(discussion_full == 1) %>%
    feols_custom(
      r2_choose_trans ~ factor(stratum_id) + r2_reliability_diff * r2_reliability_shown  + delivery_incentive_exp + factor(item_diff),
      data = .,
      fixef = c("stratum_id", "comparator_order_in_pair"),
      cluster = "group_id",
    ) %>%
    f_test("r2_reliability|item_diff", print = FALSE) %>% .$stat

  print(f2 - f1)
  return(f2 - f1)
}

set.seed(12345)
ri_fstat_diff <- conduct_ri_custom(
  test_function = compare_f_stats,
  declarations = list(ri_discussion_full),
  p = "two-tailed",
  data = r2_choices %>% filter(discuss_type %in% c("control", "discussion_full")),
  sims = 500,
  progress_bar = TRUE,
  assignments = list("discussion_full")
)

ri_fstat_diff %>% plot()

# So can I infer from this that statistical discrimination is being reduced by group discussion (rather than inducing positive taste-based discrimination?)


# Memory check ---------------------------------------------------------------

mem_check_all %>%
  bar_chart(x = discuss_type, y = mem_correct, fill = announce_after, percent = TRUE)

mem_check_all %>%
  filter(discuss_type == "choosing_only") %>%
  pull(mem_correct) %>%
  mean_na() %>%
  write_percentage("outputs/stats/mem_check_choosing_only.tex")

mem_check_all %>%
  group_by(discuss_type_label, pair_includes_trans_label) %>%
  summarise(
    mean_cl_boot(mem_correct),
    n = n()
  ) %>%
  filter(!is.na(pair_includes_trans_label)) %>%
  ggplot(aes(x = discuss_type_label, y = y, ymin = ymin, ymax = ymax, fill = pair_includes_trans_label)) +
  geom_col(position = position_dodge(0.8), width = 0.8) +
  geom_errorbar(position = position_dodge(0.8), width = 0.2) +
  scale_y_continuous(labels = scales::percent, expand = c(0, 0), limits = c(0, 1.1)) +
  theme_bw() +
  theme(    legend.position = "top") +
  labs(x = element_blank(), fill = element_blank(), y = "Mean (% correct)")



ggsave("outputs/figs/mem_check_trans_vs_non.pdf", width = 8, height = 5)

# library

mem_plot1 <- mem_check_summ %>%
  filter(outcome == "P(recalled own choices)") %>%
  ggplot(aes(x = discuss_type_label, y = estimate)) +
  geom_col(width = 0.8, position = position_dodge(0.8), fill = "#7f7f7f") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(0.8)) +
  facet_wrap(~ outcome, scales = "free_x") +
  theme_bw() +
  theme(
    plot.margin = margin(0.1, -1, 0.1, 0.1, "cm"),

  ) +
  geom_label(aes(y = 0.08, label = str_glue("{round(estimate, 2)*100}%")), position = position_dodge(0.8), size = 3, fill = "white") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1.1), expand = c(0, 0)) +
  labs(y = "Mean (95% CI)")

mem_plot1

mem_plot2 <- mem_check_summ %>%
  filter(outcome == "P(recalled other's choices)") %>%
  ggplot(aes(x = discuss_type_label, y = estimate)) +
  geom_col(width = 0.8, position = position_dodge(0.8), fill = "#7f7f7f") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(0.8)) +
  facet_wrap(~ outcome, scales = "free_x") +
  theme_bw() +
  theme(
    legend.position = c(0.5, 0.93),
    axis.text.y=element_blank(),
    plot.margin = margin(0.1, -1, 0.1, 0, "cm")
  ) +
  labs(fill = element_blank()) +
  scale_fill_manual(values = c("#FFB347", "#0096FF")) +
  geom_label(aes(group= person_category, y = 0.08, label = str_glue("{round(estimate, 2)*100}%")), position = position_dodge(0.8), size = 3, fill = "white") +

  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1.1), expand = c(0, 0))

mem_plot2

mem_plot3 <- mem_check_summ %>%
  filter(outcome == "P(recalled discussion choices)") %>%
  ggplot(aes(x = discuss_type_label, y = estimate, fill = person_category)) +
  geom_col(width = 0.8, position = position_dodge(0.8)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(0.8)) +
  facet_wrap(~ outcome, scales = "free_x") +
  theme_bw() +
  theme(
    legend.position = c(0.25, 0.93),
    axis.text.y=element_blank(),
    axis.ticks.y = element_blank(),
    plot.margin = margin(0.1, 0.1, 0.1, 0.1, "cm")
  ) +
  labs(fill = element_blank()) +
  scale_fill_manual(values = c(scales::hue_pal()(4)[[3]], scales::hue_pal()(4)[[2]], "#595959"),
                    limits = c("Discussant", "Listener")
  ) +
  geom_label(aes(y = 0.08, group = person_category, label = str_glue("{round(estimate, 2)*100}%")), position = position_dodge(0.8), size = 3, fill = "white") +

  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1.1), expand = c(0, 0))


ggarrange(
  mem_plot1 + rremove("xlab"),
  mem_plot2 + rremove("ylab") + rremove("xlab"),
  mem_plot3 + rremove("ylab") + rremove("xlab"),
  nrow = 1,
  widths = c(1.3, 1.3, 2),
  align = "v"
)

ggsave("outputs/figs/mem_check_plot.pdf", width = 8, height = 6)


mem_check_models <- list(
  "Pooled" = feols_custom(
    mem_correct ~ discuss_type * pair_includes_trans ,
    data = mem_check_all,
    cluster = "group_id",
    fixef = c("stratum_id", "video_type")
  ),
  "Split" = feols_custom(
    mem_correct ~ treat_type_r2_label + pair_includes_trans,
    data = mem_check_all,
    cluster = "group_id",
    fixef = c("stratum_id", "video_type")
  )
)


mem_check_models[[1]] %>% get_p_val("pair_includes_trans") %>%
  write_stat("outputs/stats/p_val_mem_check_diff_trans.tex", digits = 2, p_value = TRUE)

df

# Dominates / dominated [OLD] ---------------------------------------------------------------

r2_dominates %>% count_prop(pair_includes_trans, dominate_status)


r2_dominates %>%
  mutate(pair_includes_trans = factor(pair_includes_trans)) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  bar_chart(x = dominate_status, y = r2_choose_comparator, fill = pair_includes_trans, percent = TRUE, facet = discuss_type)




plot_bar_dominates <- r2_dominates %>%
  filter(!(is_listener %in% 1)) %>%
  mutate(discussion_pooled_label = ifelse(discussion_pooled == 1, "Discussion\n(pooled)", "No discussion\n(pooled)") %>%
    factor() %>% fct_relevel("No discussion\n(pooled)")) %>%
  mutate(
    pair_includes_trans_label = ifelse(pair_includes_trans, "Worker\nis trans", "Worker\nis non-trans")
  ) %>%
  bar_chart(x = dominate_status_label, y = r2_choose_comparator, percent = TRUE,
            fill = discussion_pooled_label, n_label = TRUE,
            facet = pair_includes_trans_label,
            label_size = 3.3) +
  labs(x = element_blank(), y = "Probability of selecting the worker\nin outcome round", fill = element_blank()) +
  theme_bw() +
  scale_fill_brewer(palette = "Set2", direction = -1) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = "top")


print(plot_bar_dominates)

models_dominate_status <- list(
  "Dominated" = feols_custom(
    r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(!is.na(r2_choose_trans), dominate_status == "dominated"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Neither\ndominates" = feols_custom(
    r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(!is.na(r2_choose_trans), dominate_status == "incomparable"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Dominates" = feols_custom(
    r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(!is.na(r2_choose_trans), dominate_status == "dominates"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )
)

models_dominate_status_mf <- list(
  "Dominated" = feols_custom(
    r2_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(is.na(r2_choose_trans), dominate_status == "dominated"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Neither\ndominates" = feols_custom(
    r2_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(is.na(r2_choose_trans), dominate_status == "incomparable"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Dominates" = feols_custom(
    r2_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)) %>% 
      filter(is.na(r2_choose_trans), dominate_status == "dominates"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )
)

plot_coeff_dominates <- list(
  "trans" = models_dominate_status %>% map(tidy_90) %>% bind_rows(.id = "dominate_status"),
  "non_trans" = models_dominate_status_mf %>% map(tidy_90) %>% bind_rows(.id = "dominate_status")
) %>%
  bind_rows(.id = "trans_yn") %>%
  mutate(
    pair_includes_trans_label = ifelse(trans_yn == "trans", "Worker\nis trans", "Worker\nis non-trans")
  ) %>%


  filter(str_detect(term, "discussion_pooled")) %>%
  mutate(dominate_status = fct_relevel(dominate_status, "Dominated", "Neither dominates", "Dominates")) %>%

  ggplot(aes(x = dominate_status)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "skyblue") +
  geom_errorbar(aes(ymin = conf.low_90, ymax = conf.high_90), colour = "#636363", width = 0, size = 1.2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), colour = "#636363", width = 0, size = 0.7) +
  geom_point(aes(y = estimate), show.legend = F, size = 3) +
  coord_flip() +
  facet_wrap(~ pair_includes_trans_label) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(y = str_glue("Coefficient estimate of 3-person discussion on \n probability of selecting transgender in outcome round"),
       x = element_blank())


ggarrange(
  plot_bar_dominates,
  plot_coeff_dominates,
  nrow = 2,
  align = "v",
  heights = c(2, 1),
  hjust = -5,
  vjust = 3,
  labels = "AUTO"
)

ggsave("outputs/figs/dominated_status.pdf", width = 9, height = 9)


feols_custom(
  r2_choose_comparator ~ pair_includes_trans * discussion_full * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp,
  data = r2_dominates %>% filter(discuss_type %in% c("control", "discussion_full")),
  fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id"
) %>%
  tidy_90 

# Discrimination against trans is less extreme when its dominated... (sig)
# But treatment effect of discussion is ~ invariant by dominate status (not sig)


# R1 - dominates / dominated status [OLD] ---------------------------------------------------------------

r1_dominates %>%
  mutate(pair_includes_trans = factor(pair_includes_trans)) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  bar_chart(x = dominate_status, y = r1_choose_comparator, fill = pair_includes_trans, percent = TRUE, facet = discuss_type)



plot_bar_r1_dominates <- r1_dominates %>%
  filter(!(is_listener %in% 1)) %>%
  mutate(discussion_pooled_label = ifelse(discussion_pooled == 1, "Discussion\n(pooled)", "No discussion\n(pooled)") %>%
    factor() %>% fct_relevel("No discussion\n(pooled)")) %>%
  mutate(
    pair_includes_trans_label = ifelse(pair_includes_trans, "Worker\nis trans", "Worker\nis non-trans")
  ) %>%
  bar_chart(x = dominate_status_label, y = r1_choose_comparator, percent = TRUE,
            fill = discussion_pooled_label, n_label = TRUE,
            facet = pair_includes_trans_label,
            label_size = 3.3) +
  labs(x = element_blank(), y = "Probability of selecting the worker\nin treatment round", fill = element_blank()) +
  theme_bw() +
  scale_fill_brewer(palette = "Set2", direction = -1) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = "top")

plot_bar_r1_dominates


r1_dominates %>%
  bar_chart(x = treat_type_r2_label, y = r1_choose_trans)

r1_choices %>%
  bar_chart(x = treat_type_r2_label, y = r1_choose_trans)

r2_choices %>%
  bar_chart(x = treat_type_r2_label, y = r2_choose_trans)


models_r1_dominate_status <- list(
  "Dominated" = feols_custom(
    r1_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>%   filter(!(is_listener %in% 1)) %>% 
      filter(!is.na(r1_choose_trans), dominate_status == "dominated"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Neither\ndominates" = feols_custom(
    r1_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>%  filter(!(is_listener %in% 1)) %>% 
      filter(!is.na(r1_choose_trans), dominate_status == "incomparable"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Dominates" = feols_custom(
    r1_choose_trans ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff +quality_diff,
    data = r1_dominates %>%  filter(!(is_listener %in% 1)) %>% 
      filter(!is.na(r1_choose_trans), dominate_status == "dominates"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )
)

models_r1_dominate_status_mf <- list(
  "Dominated" = feols_custom(
    r1_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(is.na(r1_choose_trans), dominate_status == "dominated"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Neither\ndominates" = feols_custom(
    r1_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(is.na(r1_choose_trans), dominate_status == "incomparable"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),
  "Dominates" = feols_custom(
    r1_choose_comparator ~ video_type_placebo + video_type_treatment + discussion_pooled + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>% filter(!(is_listener %in% 1)) %>%
      filter(is.na(r1_choose_trans), dominate_status == "dominates"),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )
)

plot_coeff_r1_dominates <- list(
  "trans" = models_r1_dominate_status %>% map(tidy_90) %>% bind_rows(.id = "dominate_status"),
  "non_trans" = models_r1_dominate_status_mf %>% map(tidy_90) %>% bind_rows(.id = "dominate_status")
) %>%
  bind_rows(.id = "trans_yn") %>%
  mutate(
    pair_includes_trans_label = ifelse(trans_yn == "trans", "Worker\nis trans", "Worker\nis non-trans")
  ) %>%


  filter(str_detect(term, "discussion")) %>%
  mutate(dominate_status = fct_relevel(dominate_status, "Dominated", "Neither dominates", "Dominates")) %>%

  ggplot(aes(x = dominate_status)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "skyblue") +
  geom_errorbar(aes(ymin = conf.low_90, ymax = conf.high_90), colour = "#636363", width = 0, size = 1.2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), colour = "#636363", width = 0, size = 0.7) +
  geom_point(aes(y = estimate), show.legend = F, size = 3) +
  coord_flip() +
  facet_wrap(~ pair_includes_trans_label) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(y = str_glue("Coefficient estimate of 3-person discussion on \n probability of selecting transgender in outcome round"),
       x = element_blank())


ggarrange(
  plot_bar_r1_dominates,
  plot_coeff_r1_dominates,
  nrow = 2,
  align = "v",
  heights = c(2, 1),
  hjust = -5,
  vjust = 3,
  labels = "AUTO"
)

ggsave("outputs/figs/dominated_status_r1.pdf", width = 9, height = 9)


feols_custom(
  r1_choose_comparator ~ pair_includes_trans * discussion_full * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp,
  data = r1_dominates %>% filter(discuss_type %in% c("control", "discussion_full")),
  fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id"
) %>%
  tidy_90 


ggarrange(
  plot_bar_dominates,
  plot_bar_r1_dominates,
  nrow = 2,
  align = "v",
  labels = c("Outcome round", "Treatment round")
)

ggsave("outputs/figs/dominated_status_combined.pdf", width = 9, height = 9)


# Dominate dominates table ---------------------------------------------------------------


models_dominates_r1_r2 <- list(
  "Outcome\nround" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * discussion_pooled * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_dominates %>% filter(!(is_listener %in% 1)),
    fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
      t = c("discussion_pooled"),
      interact = "pair_includes_trans",
      group_control = FALSE
    )
  ),
  "Treatment\nround" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans * discussion_pooled * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp + item_diff + quality_diff,
    data = r1_dominates %>% filter(!(is_listener %in% 1)),
    fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
      t = c("discussion_pooled"),
      interact = "pair_includes_trans",
      group_control = FALSE
    )
  )
)

tex_export(
  models_dominates_r1_r2,
  file = "outputs/tables/dominates_table.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id|quality_diff", "Intercept", "video", "group_control", "pair_includes_trans_alt$", control_vars, "delivery")),
  additional_header = vec_to_custom_header(
    c(" ", rep("Dep var: Chose worker (=1)", 2))
  ),
  add_rows = list_to_add_rows(
    list(
      "Controls" = c("X", "X")
    )
  )
)



# P-value on interaction terms
models_dominates_r1_r2[[1]] %>% get_p_val("pair_includes_trans:discussion_pooled:dominates") %>%
  write_stat("outputs/stats/p_val_dominates_interaction_r2.tex")

models_dominates_r1_r2[[1]] %>% get_p_val("pair_includes_trans:discussion_pooled:dominated") %>%
  write_stat("outputs/stats/p_val_dominated_interaction_r2.tex")

models_dominates_r1_r2[[2]] %>% get_p_val("pair_includes_trans:discussion_pooled:dominated") %>%
  write_stat("outputs/stats/p_val_dominated_interaction_r1.tex")

models_dominates_r1_r2[[1]] %>% get_p_val("discussion_pooled:dominated") %>%
  write_stat("outputs/stats/p_val_dominated_mf_r2.tex")


r1_dominates %>%
  filter(pair_includes_trans==1, dominated==1, discussion_pooled==1) %>%
  pull(r1_choose_comparator) %>%
  mean_na() %>%
  write_percentage("outputs/stats/prop_r1_choose_trans_dominated.tex")


r1_dominates %>%
  mutate(
    pair_includes_trans_label = ifelse(pair_includes_trans, "Worker\nis trans", "Worker\nis non-trans")
  ) %>%
  filter(phase == "phase_2", !(is_listener %in% 1)) %>%
  bar_chart(fill = pair_includes_trans_label, x = treat_type_r2_label, y = r1_choose_comparator, percent = TRUE,
            facet = dominate_status, flip = TRUE)


r1_dominates %>%
  mutate(
    pair_includes_trans_label = ifelse(pair_includes_trans, "Worker\nis trans", "Worker\nis non-trans")
  ) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  bar_chart(fill = pair_includes_trans_label, x = discuss_type, y = r1_choose_comparator, percent = TRUE,
            facet = dominate_status)


# FOcus only on costly discrimination ---------------------------------------------------------------

feols_custom(
  r2_choose_comparator ~ pair_includes_trans * discussion_full * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
  data = r2_dominates %>% filter(discuss_type %in% c("discussion_full", "control"), dominates==1, pair_includes_trans==1),
  fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id",
  lasso = TRUE,
  lasso_options = list(
    potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
    t = c("discussion_pooled"),
    interact = "pair_includes_trans",
    group_control = FALSE
  )
)

model_dominates <- feols_custom(
  r2_choose_comparator ~ pair_includes_trans * discussion_full ,
  data = r2_dominates %>% filter(discuss_type %in% c("discussion_full", "control"), dominates==1, pair_includes_trans==1),
  cluster = "group_id"
)

prop_when_dominates <- model_dominates %>% get_coeff("(Intercept)")

r2_dominates %>% count_prop(pair_includes_trans)

r2_dominates


# Split sample dominates/dominating table ----------------------------------------------------------------------------------------------------------------

r2_dominates$dominate_status

split_sample_dominates <- list(
  "Trans dominates" = feols_custom(
    r2_choose_comparator ~ discussion_full ,
    data = r2_dominates %>% filter(discuss_type %in% c("discussion_full", "control"), dominate_status == "dominates", pair_includes_trans == 1),
    fixef = c("phase", "video_type", "stratum_id", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
      t = c("discussion_pooled"),
      interact = "pair_includes_trans",
      group_control = FALSE
    )
  ),
  "Trans is dominated" = feols_custom(
    r2_choose_comparator ~ discussion_full ,
    data = r2_dominates %>% filter(discuss_type %in% c("discussion_full", "control"), dominate_status == "dominated", pair_includes_trans == 1),
    fixef = c("phase", "video_type", "stratum_id", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
      t = c("discussion_pooled"),
      interact = "pair_includes_trans",
      group_control = FALSE
    )
  ),
  "Neither dominates" = feols_custom(
    r2_choose_comparator ~ discussion_full ,
    data = r2_dominates %>% filter(discuss_type %in% c("discussion_full", "control"), dominate_status == "incomparable", pair_includes_trans == 1),
    fixef = c("phase", "video_type", "stratum_id", "comparator_order_in_pair"),
    cluster = "group_id",
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars %>% str_subset("item_diff|reliability|quality", negate = TRUE),
      t = c("discussion_full"),
      interact = "pair_includes_trans",
      group_control = FALSE
    )
  )
)

split_sample_dominates[[1]]
r2_dominates


tex_export(
  split_sample_dominates,
  file = "outputs/tables/dominates_table_split_sample.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id|quality_diff", "Intercept", "video", "group_control", "pair_includes_trans_alt$", control_vars, "delivery")),
  dep_means = list(
    "Mean: No discussion (private) + worker is non-trans" = "discuss_type == 'control' & pair_includes_trans == FALSE",
    "Mean: No discussion (private) + worker is T" = "discuss_type == 'control' & pair_includes_trans == TRUE"
  ),
  additional_header = vec_to_custom_header(
    c(" ", rep("Dep var: Chose trans worker (=1)\n(pairs with trans only)", 3))
  ),
  add_rows = list_to_add_rows(
    list(
      "Controls" = c("X", "X", "X")
    )
  )
)


model_dominates <- split_sample_dominates[[1]]



coeff_when_dominates_treat <- model_dominates %>% get_coeff("discussion_full")

(1 - prop_when_dominates) %>% write_percentage("outputs/stats/prop_when_dominates.tex", digits = 0)
(1 - prop_when_dominates - coeff_when_dominates_treat) %>% write_percentage("outputs/stats/prop_when_dominates_treat.tex", digits = 0)

coeff_when_dominates_treat %>% times_100 %>% write_stat("outputs/stats/coeff_when_dominates_treat.tex", digits = 1)
model_dominates %>% get_p_val("discussion_full") %>% write_stat("outputs/stats/p_val_only_dominates.tex", p_value = TRUE)

split_sample_dominates[[2]] %>% get_coeff("discussion_full") %>% times_100 %>% write_stat("outputs/stats/coeff_when_dominated.tex", digits = 1)
split_sample_dominates[[2]] %>% get_p_val("discussion_full") %>% write_stat("outputs/stats/p_val_only_dominated.tex", p_value = TRUE)

r2_dominates


# Effect of R1 domination on R2 ---------------------------------------------------------------

models_r1_r2_dominates <- list(
  "Just trans" = feols_custom(
    r2_choose_trans ~ (r1_n_dominates_pos + r1_n_dominated_pos) * discussion_full,
    data = r1_r2_dominates %>% filter(discuss_type %in% c("control", "discussion_full"), pair_includes_trans == 1),
    fixef = c("stratum_id", "video_type"),
    cluster = "group_id"
  ),
  "All" = feols_custom(
    r2_choose_comparator ~ (r1_n_dominates_pos + r1_n_dominated_pos) * discussion_full * pair_includes_trans,
    data = r1_r2_dominates %>% filter(discuss_type %in% c("control", "discussion_full")),
    fixef = c("stratum_id", "video_type"),
    cluster = "group_id"
  )
)


models_r1_r2_dominates[[1]] %>% tidy_90
models_r1_r2_dominates[[2]] %>% tidy_90

# r1_n_dominates:pair_includes_trans is sig, implying that (for control group as well),
# seeing a good trans option encourages you to select trans again - kind of habit formation effect / learning about "quality"

# but r1_n_dominates:pair_includes_trans:discussion_full = 0, implying that this effect is not stronger in the discussion group

# r1_n_dominated:discussion_full is sig pos >>
# implies that people that saw more dominated trans options are more likely to pick the COMPARATOR in R2, regardless of whether it's trans or not?


r1_r2_dominates %>%
  mutate(pair_includes_trans = factor(pair_includes_trans)) %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  bar_chart(x = r1_n_dominates_pos, y = r2_choose_comparator, fill = pair_includes_trans, facet = discuss_type)



# Dominates / dominating - affects group predic updating??? ---------------------------------------------------------------


group_predic_w_dominates <- group_predic_w_actual_choices %>%
  tidylog::left_join(r1_dominates_ind, by = c("ind_id", "group_id")) %>%
  filter(discuss_type %in% c("control", "discussion_full", "discussion_pair")) %>%
  mutate(discussion = discuss_type %in% c("discussion_full", "discussion_pair")) %>%
  mutate(across(c(r1_n_dominates_pos, r1_n_incomparable_pos, r1_n_dominated_pos), as.logical))


group_predic_w_dominates %>%
  bar_chart(x = r1_n_dominates, fill = discussion, y = group_predic_choose_comparator,
            facet = group_predic_includes_trans)


feols_custom(
  group_predic_choose_comparator ~ group_predic_includes_trans * discussion * (r1_n_dominates),
  data= group_predic_w_dominates %>% filter(group_predic_includes_trans),
  cluster = "group_id",
  fixef = c("stratum_id", "video_type")
) %>%
  tidy_90


group_predic_w_dominates %>%
  bar_chart(x = r1_n_dominates_pos, fill = discuss_type, y = group_predic_choose_comparator,
            facet = group_predic_includes_trans)


group_predic_w_dominates %>%
  bar_chart(x = r1_n_dominated, fill = discuss_type, y = group_predic_choose_comparator,
            facet = group_predic_includes_trans)



# No real difference in the effects discussion vs control...

# NEED A REDUCED FORM DIFFERENCE...


# CHOOSING ONLY - dominated / dominates ---------------------------------------------------------------


r1_dominates %>%
  mutate(pair_includes_trans = factor(pair_includes_trans)) %>%
  bar_chart(x = dominate_status, y = r1_choose_comparator, fill = pair_includes_trans, percent = TRUE, facet = discuss_type)


models_dominates_public <- list(
  "Treatment\nround" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans * public * (dominates + dominated) + factor(stratum_id) + delivery_incentive_exp,
    data = r1_dominates %>% filter(public==1 | control == 1) %>% filter(dominates==1, pair_includes_trans==1),
    fixef = c("phase", "video_type", "stratum_id", "delivery_incentive_exp", "comparator_order_in_pair"),
    cluster = "group_id",

  )
)

models_dominates_public[[1]] %>% get_p_val("public") %>%
  write_stat("outputs/stats/p_val_public_dominates.tex")

models_dominates_public[[1]] %>% get_coeff("public") %>%
  write_stat("outputs/stats/coeff_public_dominates.tex")


# Heterogeneity by narratives ---------------------------------------------------------------

df %>% count_prop(treat_type_r2)

feols_custom(
  as.formula(paste0(
    "r2_choose_trans ~ ",
    fml_sum(paste0(narrative_vars, "_p"))
  )),
  data = r2_with_narrative_indicators %>% filter(
    treat_type_r2 %in% c("discussion_pair", "discussion_full", "discussion_listener"),
    as.numeric(pair_includes_trans)==1
  ),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
) %>%
  tidy_90 %>%
  mutate(n = r2_with_narrative_indicators %>% filter(
    treat_type_r2 %in% c("discussion_pair", "discussion_full"),
    as.numeric(pair_includes_trans)==1
  ) %>% pull(ind_id) %>% unique() %>% length()) %>%
  mutate(narrative_val = str_replace_all(term, "(reasons_val_list_|TRUE|_p_choose|_p)", "") %>% as.numeric()) %>%
  left_join(common_narratives, by = c("narrative_val" = "reasons_val_list"), suffix = c("", "_narr")) %>%
  mutate(reasons_text_list = remove_emojis(reasons_text_list)) %>%
  mutate(
    reasons_text_label = str_glue("{reasons_text_list}\n(prob = {round(n_narr/n, 2)})") %>%
      fct_reorder(n_narr, .desc = FALSE)
  ) %>%


  ggplot(aes(x = reasons_text_label)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "skyblue") +
  geom_errorbar(aes(ymin = conf.low_90, ymax = conf.high_90), colour = "#636363", width = 0, size = 1.2) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), colour = "#636363", width = 0, size = 0.7) +
  geom_point(aes(y = estimate, colour = p.value <= 0.05), show.legend = F, size = 3) +
  coord_flip() +
  theme_minimal() +
  theme(panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(x = "Reason for choosing heard by participant in treatment round\nwhen faced with transgender choice",
       y = "Coefficient estimate on probability of\nselecting trans in outcome round")


ggsave("outputs/figs/heterogeneity_narratives.pdf", width = 7, height = 5)


# Durations ---------------------------------------------------------------


plot_durations_r1 <- choice_durations %>%
  mutate(discuss_type_label = fct_rev(discuss_type_label)) %>%
  group_by(discuss_type_label) %>%
  filter(!is.na(value)) %>%
  summarise(
    mean_cl_cluster(duration_choice_winsorised, cluster = group_id)
  ) %>%
  rename(y = estimate, ymin = conf.low, ymax = conf.high) %>%

  ggplot(aes(x = discuss_type_label, y = y, ymin  = ymin, ymax = ymax)) +
  geom_col() +
  geom_errorbar(width = 0.2) +
  coord_flip() +
  labs(
    y = element_blank(),
    x = element_blank()) +
  scale_y_continuous(limits = c(0, 75)) +
  theme_bw()

plot_durations_r1


plot_durations_r2 <- choice_durations_r2 %>%
  mutate(discuss_type_label = fct_rev(discuss_type_label)) %>%
  bar_chart(
    x = discuss_type_label, y = duration_r2_choice_winsorised, flip = TRUE
  ) +
  labs(y = "Winsorized mean duration of each choice round (seconds)",
       x = element_blank()) +
  scale_y_continuous(limits = c(0, 75)) +
  theme_bw()


ggarrange(
  plot_durations_r1,
  plot_durations_r2,
  labels = c("A. Treatment round", "B. Outcome round"),
  label.y = 0.9,
  nrow = 2,
  label.x = 0.63,
  align = "hv"
)

ggsave("outputs/figs/durations.pdf", width = 8, height = 3)


choice_durations %>%
  filter(discuss_type == "control") %>%
  pull(duration_choice_winsorised) %>%
  mean_na() %>%
  write_stat("outputs/stats/mean_duration_control.tex", digits = 1)



model_durations_r2 <- feols_custom(
  duration_r2_choice_winsorised ~ choosing_only + discussion_pooled,
  data = choice_durations_r2,
  fixef = "stratum_id",
  cluster = "group_id"
)

model_durations_r2 %>% get_p_val("choosing_onlyTRUE") %>% write_stat("outputs/stats/p_val_duration_choosing_only.tex", digits = 2)
model_durations_r2 %>% get_p_val("discussion_pooledTRUE") %>% write_stat("outputs/stats/p_val_duration_discussion.tex", digits = 2, p_value = TRUE)



mean_duration_control_r2 <- choice_durations_r2 %>%
  filter(discuss_type == "control") %>%
  pull(duration_r2_choice_winsorised) %>%
  mean_na()

model_durations_r2$coefficients["discussion_pooledTRUE"] %>% write_stat("outputs/stats/effect_duration_discussion.tex", digits = 1)
(model_durations_r2$coefficients["discussion_pooledTRUE"] / mean_duration_control_r2) %>% write_percentage("outputs/stats/effect_duration_discussion_perc.tex", digits = 0)



# Durations correlated with choices? ---------------------------------------------------------------


r2_w_durations %>% hist_basic(duration_r2_choice)

r2_w_durations$duration_r2_choice %>% quantile(c(0.25, 0.5, 0.75, 0.99))

model_durations_corr <- list(

  "Treatment round" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans * duration_choice_winsorised + item_diff + quality_diff + quality_included,
    data = r1_w_durations %>% filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase")
  ),

  "Treatment round" = feols_custom(
    r1_choose_comparator ~ pair_includes_trans * duration_choice_winsorised * discuss_type_label + item_diff + quality_diff + quality_included,
    data = r1_w_durations %>% filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase")
  ),

  "Outcome round" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * duration_r2_choice_winsorised + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_w_durations %>% filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase")
  ),

  "Outcome round" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * duration_r2_choice_winsorised * discuss_type_label + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_w_durations %>% filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase")
  )
)


tex_export(
  model_durations_corr,
  file = "outputs/tables/durations_corr.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(control_vars, c("quality_diff", "stratum_id", "Intercept", "delivery", "pair_includes_trans_alt$")),
  controls_row = TRUE,
  additional_header = vec_to_custom_header(
    c(" ", rep("Dep var: Chose worker (=1)", 4))
  )
)


model_durations_corr[[3]] %>%
  get_p_val("pair_includes_trans:duration_r2_choice_winsorised") %>%
  write_stat("outputs/stats/p_val_corr_duration_choice.tex", p_value = TRUE)


model_durations_corr[[1]] %>%
  tidy_90 



# ....... ---------------------------------------------------------------


# Heterogeneity by group characteristics ---------------------------------------------------------------

# 1. If the speakers are relatively higher on persuasive scale (within group), is it more persuasive?
# Use relative "higher" score, not absolute score, as this is random, whereas other is endogenous
r2_with_ls_listener_vs_control %>% bar_chart(x = is_listener, y = r2_choose_trans, fill = higher_ls_score_others)


# # 2. Are people who feel close to others more / less persuaded?

ls_score_var <- "higher_ls_score_others"

het_group_characs <- list(
  "No discussion (private)\n+ listeners" = feols_custom(
    as.formula(
      str_glue("r2_choose_trans ~ is_listener + {ls_score_var} + is_listener * {ls_score_var}")
    ),
    data = r2_with_relations %>%
      filter(discuss_type == "control" | is_listener == 1) %>%
      filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair")
  ),
  "No discussion (private)\n+ listeners" = feols_custom(
    r2_choose_trans ~ is_listener + close_knit_ind_others + is_listener*close_knit_ind_others,
    data = r2_with_relations %>%
      filter(discuss_type == "control" | is_listener == 1) %>%
      filter(phase == "phase_2"),
    cluster = "group_id",
    fixef = c("stratum_id", "delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair")
  ),
  "No discussion (private)\n+ 3-person discussion" = feols_custom(
    as.formula(str_glue("r2_choose_trans ~ discussion_full + {ls_score_var} +  discussion_full * {ls_score_var}")),
    data = r2_with_relations %>% filter(discuss_type %in% c("control", "discussion_full"), phase == "phase_2"),
    fixef = c("stratum_id", "delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair"),
    cluster = "group_id"
  ),
  "No discussion (private)\n+ 3-person discussion" = feols_custom(
    r2_choose_trans ~ discussion_full + close_knit_ind_others + discussion_full * close_knit_ind_others,
    data = r2_with_relations %>%
      filter(discuss_type %in% c("control", "discussion_full"), phase == "phase_2"),
    fixef = c("stratum_id", "delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair"),
    cluster = "group_id"
  )

)

het_group_characs

tex_export(
  het_group_characs,
  file = "outputs/tables/het_group_characs.tex",
  coef_rename = coef_label,
  column_widths = list(4, "1.5in"),
  coef_reorder = c("is_listener", "discussion_full", ls_score_var, str_glue("is_listener:{ls_score_var}"), str_glue("discussion_full:{ls_score_var}")),
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex(c("stratum_id", "Intercept", "video", str_subset(control_vars, "r2_reliability", negate = TRUE), "delivery")),
  dep_means = NULL,
  additional_header = c(" " = 1, vec_to_custom_header(
    rep(c("Chose trans in private outcome round\n(Phase 2 only)"), length(het_group_characs))
  )),
  stat_vec = "long",
  controls_row = TRUE
)

# Get p-value of interaction term for columns 1 and 3
het_group_characs[[1]] %>% get_p_val("is_listener:higher_ls_score_others") %>%
  write_stat("outputs/stats/pval_het_group_characs_listener.tex", digits = 2, p_value = TRUE)

het_group_characs[[3]] %>% get_p_val("discussion_full:higher_ls_score_others") %>%
  write_stat("outputs/stats/pval_het_group_characs_discussion.tex", digits = 2, p_value = TRUE)

# Get p-value for interaction term, column 4
het_group_characs[[4]] %>% get_p_val("discussion_full:close_knit_ind_others") %>%
  write_stat("outputs/stats/pval_het_group_characs_discussion_close.tex", digits = 2, p_value = TRUE)


# Look at between group heterogeneity ----------------------------------------------------------------------------------------------------------------

feols_custom(
  r2_choose_trans ~ discussion_full * high_ls_score_group,
  data = r2_with_relations %>%
    filter(discuss_type == "control" | discuss_type == "discussion_full") %>%
    filter(phase == "phase_2"),
  cluster = "group_id",
  fixef = c("stratum_id", "delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair")
) %>% list() %>% tex_export()


# Other heterogeneity ---------------------------------------------------------------


demo_vars_for_het <- demo_vars[!(demo_vars %in% c("hh_size", "hh_food_exp_pc"))] %>%
  unname() %>%
  c(., "hh_size_above_med", "hh_food_exp_pc_above_med")

het_models <- list(
  "Uninteracted term" = feols_custom(
    as.formula(str_glue(
      "r2_choose_trans ~ video_type_placebo + video_type_treatment + ({fml_sum(demo_vars_for_het)}) * discussion_full + delivery_incentive_exp + factor(team_id) + item_diff + r2_reliability_diff * r2_reliability_shown "
    )),
    data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% filter(!is.na(r2_choose_trans)) %>%
      mutate(across(c(hh_size_above_med, hh_food_exp_pc_above_med), as.numeric)),
    fixef = c("delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  ),

  "Interacted term\n(x 3-person discussion)" = feols_custom(
    as.formula(str_glue(
      "r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_full_alt + ({fml_sum(paste0(demo_vars_for_het, '_uninteract'))}) + ({fml_sum(demo_vars_for_het)}) + delivery_incentive_exp + factor(team_id) + item_diff + r2_reliability_diff * r2_reliability_shown "
    )),
    data = r2_choices_num %>%
      filter(discuss_type %in% c("control", "discussion_full")) %>%
      filter(!is.na(r2_choose_trans)) %>%
      mutate(across(c(hh_size_above_med, hh_food_exp_pc_above_med), as.numeric)) %>%
      mutate(across(all_of(demo_vars_for_het), ~.x, .names = "{.col}_uninteract")) %>%
      mutate(across(all_of(demo_vars_for_het), ~.x * discussion_full, .names = "{.col}")) %>%
      mutate(discussion_full_alt = discussion_full),
    fixef = c("delivery_incentive_exp", "video_type", "phase", "comparator_order_in_pair"),
    cluster = "group_id",
    coef_omit = "stratum_id|Intercept|video"
  )

)


tex_export(
  het_models,
  file = "outputs/tables/basic_het.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  coef_omit = vars_to_regex("\\:|stratum_id|Intercept|video|reliability|item|delivery_incentive|team_id|uninteract|discussion_full_alt"),
  controls_row = TRUE,
  stat_vec = "({std.error})",
  additional_header = vec_to_custom_header(c(" ", rep("Chose trans in\noutcome round (=1) (pairs with trans only)", 2)))
)

r2_choices_num %>% filter(
  discuss_type %in% c("control", "discussion_full")
) %>%
  pull(employer) %>% mean_na() %>%
  write_percentage("outputs/stats/perc_employer.tex", digits = 0)

het_models[[1]] %>% get_p_val("employer:discussion_full") %>%
  write_stat("outputs/stats/pval_employer.tex", digits = 2)

het_models[[1]] %>% get_p_val("employed:discussion_full") %>%
  write_stat("outputs/stats/pval_employed.tex", digits = 2)

het_models[[1]] %>% tidy_90 


het_models[[1]] %>% get_coeff("female") %>% times_100 %>% write_stat("outputs/stats/diff_mf.tex", digits = 1)
het_models[[1]] %>% get_p_val("female") %>% write_stat("outputs/stats/diff_mf_p.tex", digits = 2, p_value = TRUE)

het_models[[1]] %>% get_coeff("female:discussion_full") %>% times_100 %>% write_stat("outputs/stats/diff_mf_disc.tex", digits = 2)
het_models[[1]] %>% get_p_val("female:discussion_full") %>% write_stat("outputs/stats/diff_mf_disc_p.tex", digits = 2, p_value = TRUE)


# Phase 2 outcome round results ---------------------------------------------------------------


models_phase2_pooled <- list(
  "Chose worker in private\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans * (public_observer + discussion_pair_listener + discussion_pooled),
    data = r2_choices_num %>% filter(phase == "phase_2"),
    fixef = NULL,
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    var_types = var_types_spec,
    coef_keep = vars_to_regex(c("^pair_includes_trans$", "public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")),
    control_group = list(
      "control" = c("public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")
    )
  ),
  "Chose worker in private\noutcome round (=1)" = feols_custom(
    r2_choose_comparator ~ pair_includes_trans_alt + pair_includes_trans_alt*(video_type_placebo + video_type_treatment + public_observer + discussion_pair_listener + discussion_pooled + factor(stratum_id)) + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(phase == "phase_2"),
    fixef = c("stratum_id", "video_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    var_types = var_types_spec,
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = c("public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full"),
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    coef_keep = vars_to_regex(c("^pair_includes_trans$", "public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")),
    control_group = list(
      "control" = c("public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")
    )
  ),
  "Chose trans in private\noutcome round (=1)\n(pairs with trans only)" = feols_custom(
    r2_choose_trans ~ (pair_includes_trans) * (video_type_placebo + video_type_treatment + public_observer + discussion_pair_listener + discussion_pooled + factor(stratum_id)) + item_diff + r2_reliability_diff * r2_reliability_shown ,
    data = r2_choices_num %>% filter(phase == "phase_2") %>% filter(!is.na(r2_choose_trans)),
    fixef = c("stratum_id", "video_type", "comparator_order_in_pair"),
    cluster = "group_id",
    ri = FALSE,
    n_sims = ri_sims,
    stratum = "stratum_id",
    var_types = var_types_spec,
    lasso = TRUE,
    lasso_options = list(
      potential_controls = control_vars,
      t = c("public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full"),
      interact = "pair_includes_trans",
      group_control = TRUE
    ),
    coef_keep = vars_to_regex(c("^pair_includes_trans$", "public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")),
    control_group = list(
      "control" = c("public_observer", "public_non_observer", "discussion_pair_speaker", "discussion_pair_listener", "discussion_full")
    )
  )

)


phase_2_r2_vars_pooled <- c(
  "control_pooled",
  "public_observer",
  "discussion_pair_listener",
  "discussion_pooled"
)


p_vals_phase_2 <- models_phase2_pooled %>%
  map(
    ~ get_comparison_p_vals(.x,
                            dummy_vars = phase_2_r2_vars_pooled,
                            ri = FALSE,
                            control_group = list("control_pooled" = c("public_observer", "discussion_pair_listener", "discussion_pooled")),
                            stratum = "stratum_id", var_types = var_types_spec,
                            coef_keep = vars_to_regex(phase_2_r2_vars))
  ) %>%
  purrr::map(filter_interact_terms, vars = phase_2_r2_vars_pooled, interact = "pair_includes_trans(_alt)?") %>%
  map(filter, paste0(base, "_", term) %in% c(
    "public_observer_discussion_pair_listener",
    "public_observer_discussion_pooled",
    "discussion_pair_listener_discussion_pooled"
  )) %>%
  purrr::map(~ mutate(.x, across(c(base, term), ~ .x %>% coef_label %>% str_replace_all("\\s\\(.*\\)", ""))))

tex_export(
  models_phase2_pooled,
  file = "outputs/tables/main_phase_2_pooled.tex",
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  stat_vec = "({std.error}) [{p.value}]",
  coef_omit = vars_to_regex(c("stratum_id", "Intercept", "video", "group_control", "phase", "pair_includes_trans_alt$", control_vars, "delivery")),
  add_rows = combine_add_rows(
    list_to_add_rows(control_rows),
    add_rows = p_vals_to_add_rows(p_vals_phase_2)
  )
)


# PREF FOR WOMEN TREATMENT EFFECT ---------------------------------------------------------------

model_eff_female <- feols_custom(
  r2_choose_comparator ~ discussion_full,
  data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full"), pair_includes_female),
  cluster = "group_id",
    fixef = c("stratum_id", "video_type", "phase")
)


model_eff_female %>% get_coeff("discussion_full") %>% times_100() %>% write_stat("outputs/stats/eff_female.tex", digits = 1)
model_eff_female %>% get_p_val("discussion_full") %>% write_stat("outputs/stats/eff_female_p.tex", p_value = TRUE)

# PREF FOR WOMEN TREATMENT EFFECT - MALE PARTICIPANTS ONLY ---------------------------------------------------------------

model_eff_female_male_only <- feols_custom(
  r2_choose_comparator ~ discussion_full,
  data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full"), female == 0, pair_includes_female),
  fixef = c("stratum_id", "video_type", "phase"),
  cluster = "group_id"
)

model_eff_female_male_only %>% get_coeff("discussion_full") %>% times_100() %>% write_stat("outputs/stats/eff_female_male_only.tex", digits = 1)
model_eff_female_male_only %>% get_p_val("discussion_full") %>% write_stat("outputs/stats/eff_female_male_only_p.tex", p_value = TRUE)



# Virtue signalling - only for dominated/dominating cases ----------------------------------------------------------------------------------------------------------------

feols_custom(
  r1_choose_comparator ~ public * item_diff + public * quality_diff + discussant + item_diff + quality_diff * quality_included,
  data = r1_choices %>% filter(phase == "phase_2", !(is_listener %in% 1), !is.na(r1_choose_trans)) %>% count_prop(public, discussant),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id",
  lasso = TRUE,
  lasso_options = list(
    potential_controls = control_vars,
    t = c("public", "discussant"),
    interact = "pair_includes_trans",
    group_control = FALSE
  ),
  ri = FALSE,
  n_sims = ri_sims,
  stratum = "stratum_id",
  var_types = var_types_spec,
  coef_omit = "stratum_id|Intercept|video"
) %>%
  list() %>%
  tex_export()

r1_choices %>% filter(phase == "phase_2", !(is_listener %in% 1), !is.na(r1_choose_trans)) %>% count_prop(public, discussant) %>%
  filter(discuss_type %in% c("control", "choosing_only")) %>%
  mutate(item_diff = factor(item_diff)) %>%
  bar_chart(x = discuss_type, y = r1_choose_trans, fill = item_diff)


feols_custom(
  r1_choose_comparator ~ public * item_diff + discussant + quality_diff * quality_included,
  data = r1_choices %>% filter(phase == "phase_2", !(is_listener %in% 1), !is.na(r1_choose_trans)) %>% count_prop(public, discussant) %>%
    mutate(item_diff = factor(item_diff) %>% fct_relevel("0")) %>%
    count_prop(item_diff),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id",
  ri = FALSE,
  n_sims = ri_sims,
  stratum = "stratum_id",
  var_types = var_types_spec,
  coef_omit = "stratum_id|Intercept|video"
)


r1_choices %>% filter(phase == "phase_2", !(is_listener %in% 1), !is.na(r1_choose_trans)) %>% count_prop(public, discussant) %>%
  mutate(item_diff = factor(item_diff) %>% fct_relevel("0")) %>%
  count_prop(discuss_type, item_diff)

# only idenitfied on 84 observations...


feols_custom(
  r1_choose_comparator ~ public * dominates + discussant,
  data = r1_choices %>% filter(phase == "phase_2", !(is_listener %in% 1), !is.na(r1_choose_trans)) %>% count_prop(public, discussant) %>%
    mutate(dominates = item_diff >= 0 & quality_diff >= 0 & (quality_diff > 0 | item_diff > 0)),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair"),
  cluster = "group_id",
  ri = FALSE,
  n_sims = ri_sims,
  stratum = "stratum_id",
  var_types = var_types_spec,
  coef_omit = "stratum_id|Intercept|video"
) %>%
  list() %>%
  tex_export()


# Check heterogeneity by trans order?? ----------------------------------------------------------------------------------------------------------------

r1_trans_order <- r1_choices %>%
  filter(pair_includes_trans==1) %>%
  group_by(ind_id) %>%
  summarise(trans_sequence_r1 = paste0(round, collapse = "_")) %>%
  count_prop(trans_sequence_r1)

r2_choices %>%
  left_join(r1_trans_order, by = "ind_id") %>%
  filter(discuss_type %in% c("control", "discussion_full")) %>%
  bar_chart(x = discuss_type, y = r2_choose_trans, fill = trans_sequence_r1)

feols_custom(
  r2_choose_trans ~ video_type_placebo + video_type_treatment + discussion_full * trans_sequence_r1 + factor(stratum_id) + delivery_incentive_exp + item_diff + r2_reliability_diff * r2_reliability_shown ,
  data = r2_choices_num %>% filter(discuss_type %in% c("control", "discussion_full")) %>% filter(!is.na(r2_choose_trans)) %>%
    left_join(r1_trans_order),
  fixef = c("stratum_id", "video_type", "delivery_incentive_exp", "comparator_order_in_pair", "phase"),
  cluster = "group_id"
) %>%
  list() %>%
  tex_export()



# Are there different SOB updates depending on social closeness? ----------------------------------------------------------------------------------------------


norm_models_by_relation <- list(
  "Predicted share of people\nthat pick trans\n(community)" = feols_custom(
    n3 ~ discussion_full * relation_score_fact_z_group,
    fixef = c("stratum_id", "video_type", "phase"),
    cluster = "group_id",
    data = sn %>% 
      filter(discuss_type %in% c("control", "discussion_full")) %>%
      left_join(relations_scores_group, by = "group_id") %>% 
      ungroup %>% 
      mutate(relation_score_fact_z_group = as.numeric(relation_score_fact_z_group > median_na(relation_score_fact_z_group))),
    ri = FALSE
  ),
  "Predicts that other picks trans (=1)\n(within group-of-3)" = feols_custom(
    fml = group_predic_choose_trans ~ discussion_full * relation_score_fact_z_group + item_diff + reliability_diff + reliability_shown,
    fixef = c("stratum_id", "video_type", "phase"),
    cluster = "group_id",
    data = group_predic %>% 
      filter(discuss_type %in% c("control", "discussion_full")) %>%
      filter(!is.na(group_predic_choose_trans)) %>%
      left_join(relations_scores_group, by = "group_id") %>% 
      ungroup %>% 
      mutate(relation_score_fact_z_group = as.numeric(relation_score_fact_z_group > median_na(relation_score_fact_z_group))),
    ri = FALSE
  )
)

tex_export(
  norm_models_by_relation,
  coef_omit = vars_to_regex("stratum_id|Intercept|video_type|n1|n2|b11|item_diff|reliability_", control_vars),
  coef_rename = coef_label,
  gof_map = fe_label_no_fe,
  dep_means = list("Mean: No discussion (private)" = "discuss_type == 'control'")
)

norm_models_by_relation[[2]] %>% get_p_val("discussion_full:relation_score_fact_z_group") %>% write_p_val("outputs/stats/norms_by_relation_p.tex")



# GRoup predic accuracy by relation ----------------------------------------------------------------------------------------------

group_predic_w_actual_choices %>%
  filter(discuss_type %in% c("control")) %>%
  left_join(relations_scores_group, by = "group_id") %>% 
  filter(group_predic_includes_trans) %>% 
  feols_custom(correct_predic ~ relation_score_fact_z_group, data = ., cluster = "group_id", fixef = "stratum_id") %>% 
  get_p_val("relation_score_fact_z_group") %>% 
  write_p_val("outputs/stats/group_predic_accuracy_by_relation_p.tex")
